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SUMMARY 

A method is proposed for accurately describing arbitrary-shaped free boundaries in finite- 
difference schemes for elastodynamics, in a time-domain velocity-stress framework. The ba- 
sic idea is as follows: fictitious values of the solution are built in vacuum, and injected into 
the numerical integration scheme near boundaries. The most original feature of this method 
is the way in which these fictitious values are calculated. They are based on boundary condi- 
tions and compatibility conditions satisfied by the successive spatial derivatives of the solution, 
up to a given order that depends on the spatial accuracy of the integration scheme adopted. 
Since the work is mostly done during the preprocessing step, the extra computational cost is 
negligible. Stress-free conditions can be designed at any arbitrary order without any numeri- 
cal instability, as numerically checked. Using 10 grid nodes per minimal S-wavelength with a 
propagation distance of 50 wavelengths yields highly accurate results. With 5 grid nodes per 
minimal S-wavelength, the solution is less accurate but still acceptable. A subcell resolution of 
the boundary inside the Cartesian meshing is obtained, and the spurious diffractions induced 
by staircase descriptions of boundaries are avoided. Contrary to what occurs with the vacuum 
method, the quality of the numerical solution obtained with this method is almost independent 
of the angle between the free boundary and the Cartesian meshing. 

Key words: Free surface. Seismic modeling. Velocity-stress formulation. Numerical methods. 
Finite-difference methods, ADER schemes. Boundary conditions. Compatibility conditions. 



1 INTRODUCTION 

Various approaches have been proposed for simulating the prop- 
agation of elastic waves with free boundaries. The first app roach 
is ba sed on variational methods , as done in finite elements jPavl 
Il977ii) . spectral finite elements (Komatitsch & Vilottel [1998') and 
discontinuous Galerkin (Ben Jemaa et al. , 2 007|) . These methods 
provide a fine geometrical description of boundaries by adapting 
the mesh to the boundaries. Boundary conditions are accounted for 
weakly by the underlying variational formulation. However, a grid- 
generating tool is required, and small time steps may result from 
the smallest geometrical elements and from the stability condi- 
tion. T he SAT methods based on energy estimates jCarpenter et all 
1 1994 ) avoid these limitations by introducing Cartesian grids and 
give time-stable high-order schemes with interfaces. However and 
up to our knowledge, these methods have not been applied so far to 
elastodynamics with free boundaries. 

The second approach used in this context is based on the 
strong form of el astodynamics, as done in finite differences and 
spectral methods jTessmer & Koslo3 . ll994l) . In seismology, finite 
differences are usually implemented on staggered Cartesian grids, 
either with completely staggered stencils (CSS) or with the recently 



developed partially staggered stencils (PSS). With CSS, the veloc- 
ity and st ress componen ts are distributed between different node 
positions jVirieuxl 11983) . With PSS, all the velocity components 
are computed at a single node, as are the stress components, al- 
though the latter are shifted by hal f a node in two separate g rids. 
Second-order ( Saenger et al., 2000"; Saenger & Bohlen, 20 (341) and 
fourth -order tBohlen & Saenger . .2003; Cruz-Atienza & Virieuxl 
|2004) spatially-accurate PSS have been developed; for further dis- 
cussion, we denote them PSS-2 and PSS-4, respectively. Unlike 
variational methods, finite differences require special care to incor- 
porate the free boundary conditions strongly. There exist two main 
strategies for this purpose: 



(i) First, the boundaries can be taken into account implicitly 
by adju s ting t h e physical parame ters locally iKellv etaU 119761: 
IVirieuxl 1 19861 : iMuir et all Il992h . The best-known implicit ap- 
proach is the so-called vacuum method (Za hradn^ . 1 1 9951 : iGravesl 
1996; Moczo et al., 2002; Gehs et al., 2005). For instance, the vac- 
uum method applied to PSS involves setting the elastic parameters 
in the vacuum to zero, and using a small density value in the first ve- 
locity node in the vacuum to avoid a division by zero. However, this 
easy-to-implement method gives at best second-order spatial accu- 
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racy. In addition, a systematic numerical study has shown that the 
accuracy of the solution decreases dramatically when the angle be- 
tween the boundary and the meshing increases jBohlen & SaengerL 
L2OO6). Lastly, applying the vacuum method sometimes gives ri se to 
instabilities: see for instance PSS-4 teohlen & Saengei],l2003h . 

(ii) A secon d idea is to expli citly change the scheme near 
the boundaries jKellv et all 1 19761) . The best-known explicit ap- 
proach is the so-called image method, which was deve loped for 
dealing with flat boundaries to fourth-order accu racv tLevandei . 
19881) and then ex t ended to variable topo graphies jjih et al.Lll98j 



O 



Robertssonl. 1 1994 |Zhang & Cheij I2OO6I) . However image meth- 



ods require a fine grid to reduce the spurious diffractions up to an 
acceptable level. To avoid this spatial oversampling, various tech- 
niques have bee n proposed, such as grid refinement in the vicinity 
of the boundary ( lRodrigueslll995h or a djusted finite-difference ap- 
proximations: see ( lMoczoetal.Ll2007h for a review on these sub- 
jects. 



The aim of this paper is to present a finite-difference approach 
accounting for free boundaries without introducing the aforemen- 
tioned drawbacks of the vacuum and image methods. The require- 
ments are as follows: smooth arbitrary-shaped boundaries must be 
treated as simply as straight boundaries; the accuracy of the method 
must not depend on the position of the boundary relative to the 
meshing; and lastly, the computations must be stable even with 
very long integration times. We establish that these requirements 
can be met by applying an explicit approach involving fictitious 
values of the solution in the vacuum. In previous studies, inter- 
face problems in t he context of elasto d ynam i cs were investigated 
in a similar way (^ Piraux & Lombard 1200 ll ; [Lombard & Pirau3. 



12004. 2006). The fictitious values are high-order Taylor expansions 
of the boundary values of the solution. Estimating these bound- 
ary values involves some mathematical background, in order to 
compute the high-order boundary conditions; to determine a mini- 
mal set of independent boundary values; lastly, to perform a least- 
square numerical estimation of this minimal set. To help the reader, 
subroutines in FORTRAN are proposed freely at the web page 
http : //w31ma. cnrs-mrs . f r/~ Ml/Software/. These subrou- 
tines enable a straightforward implementation of the algorithms de- 
tailed in the present paper. 

The disadvantage here is that the above requirements can- 
not be fully satisfied if staggered-grid schemes are used. Single- 
grid finite-difference schemes are therefore chosen, where all the 
unknowns are computed at the same grid nodes. Our numerical 
experiments are based on the hi gh-order ADER schemes w hich 
are widely used in aeroacoustics jSchwartzkopff et^lll2005l) . Al- 
though t hese schemes are not yet widely used in the field of seis- 
mology jPumbser & Kaseii 12006^. they have also great qualities 
because of their accuracy and their stability properties: using 10 
grid nodes per minimal S-wavelength with a propagation distance 
of 50 wavelengths gives highly accurate results. Moreover, on 
Cartesian grids, these methods do not require much more compu- 
tational memory than staggered-grid schemes. 

This paper is organized as follows. Section 2 deals with the 
continuous problem: the high-order boundary conditions and com- 
patibility conditions are stated. These conditions are useful for han- 
dling the discrete problem presented in section 3, where the focus is 
on obtaining fictitious values of the solution in the vacuum. In sec- 
tion 4, numerical experiments confirm the efficiency of this method 
in the case of various topographies. In section 5, conclusions are 
drawn and some prospects suggested. 



Figure 1. Boundary F between a solid and vacuum. 

2 THE CONTINUOUS PROBLEM 
2.1 Framework 

Let us consider a solid fl separated from the vacuum by a boundary 
r (Figure [Til. The configuration is in-plane and two-dimensional, 
with a horizontal axis x and a vertical axis z pointing respectively 
rightward and downward. The boundary F is described by a para- 
metric expression {x{t), z{t)) where the parameter r describes 
the sampling of the boundary. The tangential vector and the nor- 
mal vector are t — '^{x (r), z (r)) and n — "^{—z (r), x (r)), 
respectively, with x (r) = ^(r), z (r) = 77 (t), and T refers 
to the transposed vector. We assume the spatial derivatives at any 
point of the boundary to be available, as specified below. 

The solid Q is assumed to be linearly elastic, isotropic, and to 
have the following constant physical parameters: the density p and 
the Lame coefficients A, ft,. The P- and S-wave velocities are Cp — 
J p./ p. A velocity-stress formulation is 



^(A + 2p)/ p and Cs 
adopted, hence the unknowns are the horizontal and vertical com- 
ponents of the elastic velocity d = '^{vx, Uz), and the independent 
components of the elastic stress tensor cr — {a^x, o^z, o-zz)- 
Setting 
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2.2 High-order boundary conditions 

At any point P(r) on the free surface F (Figure[T), the stress tensor 
satisfies the homogeneous Dirichlet conditions cr.n — 0. These 
zero-th order boundary conditions are written compactly 

L°{r)U"{x{r),zir),t) = 0, (2) 

where U^' is the limit value ofU at P and L° is the matrix 

-z'(t) x'{t) 
-z'{t) x'{t) 

From now on, the dependence on r is generally omitted. To de- 
termine the boundary conditions satisfied by the first-order spatial 
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derivatives of U, two tasks are performed. First, the zeroth-order 
boundary conditions ^ are differentiated in terms of t. The time 
derivative is replaced by spatial derivatives using the conservation 
laws ([TJ, which gives 



\ ox oz J 



(3) 



Secondly, the zeroth-order boundary conditions l[2j are differenti- 
ated in terms of the parameter r describing V. The chain-rule gives 

Vdr / V dx dz ) 

Since the matrix dV-^ /dr in ^ involves x and z , it accounts 
for the curvature of F at P. Setting the block matrix 








equations l|2j, l[3j and l|4j give the boundary conditions up to the 
first-order 



f = 0, 



with 



U ^ hm [ U,— U, U 

A/£n-.p V ox oz 

Let fc > 1 be an integer whose value will be discussed in section 
[3l To get the boundary conditions up to the fc-th order, one deduces 
from ([2J 



Qk 



a — 0, k. 



(5) 



The r-derivatives are replaced by spatial derivatives by applying 
{k — a)-times the chain rule. The t-derivatives are replaced by spa- 
tial derivatives by injecting Q-times the conservation laws (TJ. The 
boundary conditions so-obtained up to the k-th order can be written 
compactly 



with 



0, 



lim 



dx° 



'dz^ dz" 



(6) 



,(7) 



where a = 0, k and /3 = 0, a. The vector U has n„ = 
5(A: + 1) (A: + 2)/2 components. L'' is a n; x n„ matrix, with 
ni = (k+l) (fc+2). This matrix involves the successive derivatives 
of the curvature of F at P. Computing L'^ with fc > 2 is a tedious 
task, which can be greatly simplified by using computer algebra 
tools. 



2.3 Compatibility conditions 

The second spatial derivatives of stress components are linked to- 
ether by th e compatibility condition of Barre-de Saint Venant 
Lovelll944h 



dxd z 
with 



■ a2 ■ 



d (Tx 



+ a-i ■ 



d (Tx 



+ 0L2 



(8) 



A + 2/i 
4(A + m)' 
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This compatibility condition is a necessary and sufficient condition 
for the strain tensor to be symmetrical. If fe > 2, it can be dif- 
ferentiated (fc — 2)-times in terms of x and z. With fc > 2, one 
obtains — k (k — l)/2 relations; with fc < 2, nc = 0. Unlike 
the boundary conditions, these compatibility conditions are satis- 
fied everywhere in Q,: in particular, they are satisfied at P on F. The 

vector of boundary values can therefore be expressed in terms 

- fc 

of a shorter vector U with n^, — rtc independent components 

A n algorithm for building th e n„ x (n„ — ric) matrix G'' is given 
in lLombard & Pirauxl ( [2006h . 



3 THE DISCRETE PROBLEM 
3.1 Numerical sclieme 

To integrate the hyperbolic system (TJ, we introduce a single Carte- 
sian lattice of grid points: (xi,Zj,tn) = {ih,jh, 71 A t), where /i 
is the mesh spacing and A f is the time step. Unlike with staggered 
grids, all the unknowns are computed at the same grid nodes. The 
approximation U"j of U{xi, Zj ,tn) is computed using any ex- 
plicit, two-step, and spatially-centred finite-difference scheme. A 
re view of the huge bo dy of literature on fin ite-differences is given 
in lLeVeaud h992h and lMoczo etai] ilOOj) . 

Here we propose to use ADER schemes, that allow to 
reach easily arbitrar y high -order of time and space accuracy 
jSchwartzkopif et all 12005 ). On Cartesian grids, these finite- 
volume integration schemes originally developed for aeroacous- 
tic applications are equival ent to finite-difference L ax-Wendroff- 
type integration schemes l lLorcher & Munzl . |2005|) . In the nu- 
merical experiments described in section |4l we use a fourth- 
order ADER integration scheme. This scheme is stable under the 
Courant-Friedrichs-Lewy (CFL) condition CpAt/h < 0.9 in 
2D; as usually with single- grid schemes, it is slightly dissipative 
dSchwartzkopff et al.L[20o3) . 

Many other single-grid schemes can be used in this context. 
In particular, the method described in the next sub sections has been 
successfully combined with flux-limiter schemes ( lLeVeauj . ll992l) 
and with the standard second-order Lax-Wendroff scheme. Diffi- 
culties have been encountered with dissipative-free schemes based 
on centred staggered-grid finite-difference schemes, as we will see 
in section[331 



3.2 Use of fictitious values 

Time-marching at grid-points where the stencil crosses F requires 
fictitious values of the solution in the vacuum, which have to be 
determined. The question arises as to how to compute, for instance, 
the fictitious value U} j at the grid point {xi, zj) in the vacuum, 
as sketched in Figure (2] Let P(r) be the orthogonal projection of 
{xi, zj) on F, with coordinates {xp — x{t), zp — z{t)). At any 
grid point {xi, Zj), we denote 



/5, 



{x,~xpr-''{ 



zp 



zp 



fc! 



the 5 X matrix containing the coefficients of fc-th order Taylor 
expansions in space at P, where I5 is the 5x5 identity matrix, 
a = 0, fc, and {3 — 
as the Taylor-like extrapolation 



a. The fictitious value U*j j is defined 
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Figure 2. Determination of the fictitious value U*jj required for time- 
marching at neighboring grid nodes. P is the orthogonal projection of 
{xj, zj) on r. The rip grid nodes in Q and inside the circle C centred 
at P with a radius d are denoted by +. 

U'i^j = u]^jU\ (10) 
where defined by Q still remains to be estimated. 

3.3 Reduced vector of boundary values 

Before determining in dlOt . we first reduce the number of inde- 
pendent components it contains. The expressions obtained in sec- 
tion |2] are used for this purpose. The linear homogeneous system 
following from ([6} and ([9} is 

L'=G*C/*=0. (11) 

This system has fewer equations {ni) than unknowns (n„ — ric). It 
therefore has an infinite number of possible solutions that constitute 
a space with the dimension riv ~ ~ ni. Let K i^kQk be a {jiv — 
ric) X (n„ — He — ni) matrix containing the basis vectors of the 
kernel of G'' . The general solution of dl lb is therefore 

U" =K^ka''U\ (12) 

where the Uv — ric ~ rii components of [7* are real numbers. In- 
jecting l ll2t into ^ gives 

U'' = G''K^kGkU''. (13) 

The computation of K i^kQk is a key point. For this purpose, we 
use a classical linear algebra tool: singular value decomposition 
of G^. Technical det ails can be found in the Appendix A of 
[Lombard & Pirauxl ( l2004h . 

3.4 Computation of fictitious values 

Let us now consider the rip grid points of Q in the circle C centred 
at P with a radius d; for instance, Tip = 8 in Figure |2] At these 
points, we write the fc-th order Taylor expansion in space of the 
solution at P, and then we use the expression l ll3t . This gives 

_ (14) 

= nl^G'^K^kakU'' +0{h!=+'-). 

The set of rip equations l ll4b is written compactly via a 5 Wp x (uv — 
ric — ni) matrix M 



(?7(., t„))c =M[7' + 0(fe"+i), (15) 

where {U{., tn))(< is the vector containing the exact values of the 
solution at the grid nodes of Q, inside C. These exact values are re- 
placed by the known numerical values {U")^, and Taylor rests are 
removed. From now on, numerical values and exact values of the 
fields are used indiscriminately. The discrete system thus obtained 
is overdetermined (see the remark (i) about d and typical values of 
Hp in subsection l3.5l >. We now compute its least-squares solution 

U'' = M-^ {U")c , (16) 

where the (n„ — — n;) x 5 rip matrix denotes the pseudo- 

inverse of M. From dlOt . l ll3t and l ll6t . the fictitious value in the 
vacuum at (xi, zj) is 

u},j = nt/G'-Fs:i.G.M-i([7")c 

(17) 

= A... ([/"),. 

The 5 X 5np matrix is called the extrapolator at (xj, zj). 
The fictitious values have no clear physical meaning. They only 
allow, by interpolation with numerical values inside fl, to recover 
the high-order Dirichlet conditions l|7j. 

3.5 Comments and practical details 

The extrapolation method described in section [T4l has to be applied 
at each grid point (/, J) in the vacuum where a fictitious value 
is required for the time-marching procedure. Useful comments are 
proposed about this method: 

(i) The radius d of C must ensure that the number of equations 
in ( list is greater than the number of unknowns: 

e{k, d) = > 1. (18) 

riv — ric — rii 

No theoretical results are available about the optimal value of e. 
However, numerical studies have shown that a definite overestima- 
tion ensures long-term stability: typically, £ ~ 4. Various strategies 
can be used to ensure i fTSt , such as an adaptative choice of d de- 
pending on the local geometry of F at P. Here we adopt a simpler 
strategy consisting in using a constant radius d. With fc = 3, nu- 
merical experiments have shown that d = 3.2 /i is a good candidate 
for this purpose. In this case, one typically obtains rip « 15. 

(ii) Since the boundary conditions do not vary with time, the 
extrapolators A/,,/ in l ll7b can be computed and stored during a 
pre-processing step. At each time step, only small matrix-vector 
products are required. 

(iii) The extrapolators A/,j account for the local geometry of 
r at the projection points P on F via (section |2^ . Moreover, 
they incorporate the position of P relative to the Cartesian mesh- 
ing, via Ilij J14t and 11/, j The set of extrapolators therefore 
provides a subcell resolution of F in the meshing, avoiding the spu- 
rious diffractions induced by a naive description of the boundaries. 

(iv) The stability of the method has not been proved. However, 
numerical experiments clearly indicate that the CFL condition of 
stability is not modified compared with the case of a homogeneous 
medium. The solution does not grow with time, even in the case of 
long-time simulations (see section l43] ). 

(v) In a previous one-dimensional study jPiraux & Lombarcl 
l200lh . the local truncation error of the method has been rigorously 
analysed, leading to the following result: using the fictitious values 
dlVI l ensures a local r-th order spatial accuracy if fc > r, where 
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Figure 3. Staggered-grid schemes with a plane boundary T parallel to the 
meshing: two cases can be distinguished, depending on the position of F 
relative to the meshing. Case (a), where the fictitious stress is estimated, 
works well, while case (b), where the fictitious velocity is estimated, leads 
to long-term instabilities. 



Q 



Similar behavior is observed with PSS-4, but after a longer time: 
the numerical solution generally works well during a few thousand 
time steps, before growing in a unstable manner. 

The extrapolation method presented here is therefore not rec- 
ommended for use with staggered-grid schemes, especially PSS-2, 
except in the trivial case sketched in Figure[3]-(a). 



r is the order of spatial accur acy of the scheme. In 2D confi gra- 
tions with material interfaces ( [Lombard & Pirau^l2004l[2006l) . no 
proof has been conducted, but numerical experiments have shown 
that the r-th order overall accuracy is also maintained by taking 
k = r. Note that a slightly smaller order of extrapolation can be 
used: k = r — \ suffices to provide r-th order overall accuracy 
( lGustafssorill975l) . The value fc = 3 is therefore used for the 
fourth-order ADER scheme. 

(vi) The extrapolators do not depend on the numerical scheme 
adopted. They depend only on k and on physical and geomet- 
rical features. Standard subroutines for computing the extrapo- 
lators A], J can therefore be developed and adapted to a wide 
range of schemes. Subroutines of this kind are freely available 
in FORTRAN at the web page |http: //w31nia. cnrs-mrs .fr/'l 
Ml/Software/. 

3.6 Case of staggered-grid schemes 

Instead of using a single-grid scheme as proposed in section [3T1 
readers may be interested in adapting our approach to staggered- 
grid schemes such as CSS or PSS (see section[T]for the definition 
of these terms). However, in the case of some of the boundary po- 
sitions relative to the meshing, computational instabilities occur, 
especially when long-time integration is considered. 

To understand why this is so, let us consider PSS-2. Taking 
a simple flat boundary to exist between the medium and the vac- 
uum leads to two typical geometrical configurations. At one po- 
sition of the free surface, the boundary discretization will require 
only the stress field to be extrapolated (Figure[3]-(a)). Our procedure 
works satisfactorily with this type of discretization at any order k. 
It also yields stable and accurate solutions when dealing with PSS- 
4, contrary to the vacuum method. Using 10 grid nodes per mini- 
mal S-wavelength gives similar performance in this case to those 
of our numerical experiments based on the ADER scheme, which 
are shown in section|4] 

At another position of the free surface where only extrapo- 
lated velocities are required within a wide zone (Figure[3]-(b)), our 
procedure results in instabilities. The reason for this problem is as 
follows: fictitious velocities involve first-order boundary conditions 
([3} and higher-order conditions (see section lZ2l (. but they do not in- 
volve the fundamental zeroth-order Dirichlet conditions ([2}. Since 
the latter conditions are never enforced, an increasing oscillating 
drift occurs near the boundary, which invalidates the computations. 



3.7 Case of non-smooth geometries 

Up to know, we have assumed that the boundary T was sufficiently 
smooth at the projection points, being at least C^'^^ at each P, 
where fc > is the order of differentiation defined in section 
13.51 Let us assume now that F is only at a point P, with 
K < fc -I- 1. Then, the components of L'^ in ^ involving the 
derivatives ^(t) and '^{r) (q = if + 1, fc + 1) of the 
parametric representation are discontinuous, invalidating locally 
the method proposed. In our software, we have implemented the 
following rough treatment: 

(i) \f K = Q, the boundary owns a corner and the solution has 
an integrable singularity. The comer is replaced by an arc of circle 
centred at Q with radius 5 (figure|4j, leading to a curve. 

(ii) If < A' < fc + 1, as in the previous case at Pq and Pi with 
K = I, the values of ^(r) and ^(r) (q = A' + 1, fc + 
1) are taken indiscriminately on one side or the other of the point 
considered. 

No numerical instabilities were observed if 5 (in case (i)) or the 
radius of curvature (in case (ii)) is much greater than h. It is agreed 
that the accuracy of computations is no more controlled in the cases 
(i) and (ii), especially the convergence towards the exact solution. 

More sophisticated treatments of geometrical singularities, 
such as space-time mesh refinement (Berger & LeVeque, 1998), re- 
quire further investigation, which is out of the scope of the present 
paper New studies are also needed in the case of merging bound- 
aries, occuring for inst ance when an inter nal material interface 
reaches the free surface jMoczo et aLll2004l) . 



4 NUMERICAL EXPERIMENTS 
4.1 Configurations 

The time evolution of the source is a Ricker wavelet 

g{t) = {2{nU{t~ U)f - 1} e- ^= - , (19) 

where fc is the central frequency, and tc = 1/ fc- The maximal 
frequency /max defined by \g{fuia.^) /g{fc)\ = 0.5 (the tilde des- 
ignates the Fourier transform) is /max ~ 1-6 fc- We will adopt this 
frequency /max for our rule of thumb about the number of grid 
nodes per S-wavelength. The following values of the physical pa- 
rameters will be used in all the following tests: p = 2400 kg/m^. 



6 B. Lombard, J. Piraux, C. Gelis, J. Virieux 

(a) 

14 
12 
10 

8 

6 

4 



^ %, W 1» 3St. 

(b) 





(a) 



-5E-5 



zooml 


[ 




zoom2 


zoomS 




r 






0.5 


1.5 


2 2.5 






t(s) 








(b) 





1.5E-4- 




-I nF.-4- l ^ ^ ^ ^ 

0.3 0.4 0.5 0.6 0.7 

t(s) 



(C) 




1.05 1.2 1.35 1.5 1.65 

l(s) 



(d) 




1.8 1.95 2.1 2.25 24 

t(s) 



Figure 5. Test 1: snapshots of Vx at the initial instant (a), at mid-term (b) 
and at the final instant (c). 



Figure 6. Test 1: time history ofvx (a). Zooms on successive time windows, 
with various discretizations (b,c,d): the number after # denotes the number 
of grid nodes per minimal S-wavelength. 
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Figure 8. Test 2: parametric study of the relative error in terms of the 
boundary 's angle 6, with various discretizations. The number after # de- 
notes the number of grid nodes per minimal S-wavelength. 



t(s) 



Figure 7. Test 2: snapshot ofvz at the final instant (a). Numerical and exact 
time histories of (b). 



Cp = 4500 m/s, and Cs — 2200 m/s. Lastly, the mesh size and the 
time step satisfy Cp At/h — 0.85. 

The simulations are performed on a PC Pentium at 3 GHz 
with 2 GB of RAM. The results of tests 1 and 2, with constant 
and null curvature of F, compare quantitatively with analytical so- 
lutions denoted by a solid line. Test 3, with a variable curvature, is 
purely qualitative. Test 4 shows the slow decrease in the mechani- 
cal energy which occurs during very long integration times, which 
confirms the stability of the method. 



4.2 Test 1: circular boundary 

Computational efficiency. Let us consider a circular cavity contain- 
ing vacuum, with radius 1 km, at the center of a 18 km x 18 km 



domain. In a first part, the mesh spacing is ft = 25 m. The source 
is a rightward-moving plane wave, with /max = 8 Hz, ensuring 22 
grid nodes per minimal P- wavelength and 10 grid nodes per mini- 
mal S-wavelength at that frequency. 

During the pre-processing step, the program finds the 616 grid 
nodes where fictitious values are required; it also computes and 
stores the 616 extrapolators defined by the expression (TTJ. Time 
integration is performed in 550 time steps, which corresponds to a 
propagation time of 2.75 s and a propagation distance of 22 mini- 
mal wavelengths. The preprocessing step takes 21 s of CPU time. 
The time integration takes 1 100 s of CPU time, including 28.6 s in- 
duced by the computation and by the use of fictitious values, which 
amounts to an extra time cost of only 2.6 %. Figure[5] shows snap- 
shots of Vx at the initial instant (a), after 275 time steps (b) and after 
550 time steps (c). 

Quantitative study. In a second part, three discretizations are 
considered: h — 25 m, 50 m, 60 m, corresponding respectively to 
10, 5 and 4 grid nodes per minimal S-wavelength. A receiver is set 
just above the cavity, at the position (9 km, 10.2 km), that can be 
seen on Figure [5j it mainly records the waves propagating along 
the boundary. Numerically, these waves are highly sensitive to the 
quality of the fictitious values defined in section [3] 

Figure|6l-(a) shows the time history of Vx at the receiver. In this 
time window, three main wave packets are generated; with the scale 
of Figure [6]-(a), the third packet cannot be seen. The amplitude is 
divided by a factor of approximately 30 from one packet to the fol- 
lowing one. For the sake of clarity, zooms around each wave packet 
are shown in Figure[6l-(b,c,d). These solutions are compared with 
an exact solution computed thanks to inverse Fourier transforms on 
4096 frequencies, with 1.25 10~^ Hz as the sampling frequency; 
each harmonic component is expanded into 60 Bessel modes. The 
agreement between the numerical and the analytical values is very 
good when 10 grid nodes per wavelength are used, even at very 
small amplitudes (d). For 5 and 4 grid nodes per wavelength, the 
solution is slightly less accurate, but it is still acceptable. 
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Figure 9. Test 3: snapshots ofvz, for various sinusoidal topographies. 



4.3 Test 2: inclined straight boundary 

The Garvin's problem. As a second test, we take a plane boundary 
inclined against the Cartesian mesh. The domain under investiga- 
tion is 18 km wide and 12 km high, with the origin of the coor- 
dinates on the top and left. The mesh spacing is h — 10 m. Four 
receivers at (10 km, 1.8 km), (11 km, 1.6 km), (12 km, 1.4 km) and 
(13 km, 1.2 km) belong to the free boundary which is inclined at an 
angle of 6* — 11.3 relative to Ox. 

An explosive source S is buried at (9 km, 2.1 km), with 
/max = 24 Hz. The distance between the source and the free sur- 
face is roughly 100 m < Ap/3, where Ap is the wavelength of the 
compressional waves at frequency /c, and hence large Rayleigh 
waves are generated, with a velocity Cr — 2054 m/s. To prevent 
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Figure 10. Test 4: snapshot of Vz (a) and time history of the mechanical 
energy (b). 



the occurence of spurious oscillations, the source is spread out nu- 
merically over a radius of i? = Ap/7.5 = 40 m. The source is 
weighted by a gaussian law with a standard deviation R/2 — 20 
m. The spatial discretization ensures the sampling of roughly 18 
grid nodes per minimal P-wavelength, and 9 grid nodes per mini- 
mal S-wavelength at the frequency /max- 

Figure |7]-(a) shows a snapshot of Vz after 1200 time steps, 
corresponding to a propagation time of 2.25 s and a propa- 
gation distance of 55 minimal wavelengths. Direct cylindrical 
waves are observed, together with converted PP waves, con- 
verted PS waves (with an almost linear wavefront), and Rayleigh 
waves. In Figure |7]-(b), the time history of Vz recorded at the re- 
ceivers can be favourably compared with an exact solution. The 
latter is obtained by convolving the Green's fu nction o b tained 
by the well-known Cagniard-de Hoop method i lGarvinl 1 19551 : 
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ISanchez-Sesma & Iturraran-Viverosl . l2006l) with the source wavelet 
il9\ and with the discrete source spreading. 

Influence of the slope. To quantify the effects of the angle be- 
tween the boundary and the Cartesian meshing on the numerical 
solution, we perform a parametric study of the error in terms of 6. 
Ten angles are considered, from 9 = 0° to 6 = 45° in steps of 5°. 
In each configuration, the waves are measured at the free boundary 
after propagating for 65 minimal wavelengths. The error of v.n is 
measured in norm L2, and then it is normalized by the norm L2 of 
the exact time history of w.n. 

The results of this study are shown in Figure [8] with vari- 
ous discretizations: 5, 10, 20 grid nodes per minimal S-wavelength. 
With a given h, the error is almost constant and independent of 6. 
This constitutes a crucial advantage of our method over the vacuum 
method, where the error at 45° is much greater than that at 0° : it 
means that an extremely fine discretization is required to obtain 
accurate results with the v acuum method when arbi trary- shaped 
boundaries are encountered jBohlen & Saengeill2006h . 

4.4 Test 3: sinusoidal boundary 

Since boundaries not related to the finite-difference grid can be in- 
cluded, the third test is performed on a sinusoidal free boundary, 
with a peak-to-peak amplitude of 800 m and various wavelengths: 
0.5 km, 1 km and 2 km. The sinusoidal curve is centred around 
2 = 1 km. The source S is located at (9 km, 1 km). The other pa- 
rameters are the same as in test 2. Figure |9] shows snapshots of 
at the final instant 2.25 s. One can clearly see how the wavelength 
of the sinusoidal boundary influences the diffracted fields. 

Convergence studies (not shown here) were performed in 
these three cases, by comparing solutions computed on finer grids. 
We again concluded that accurate solutions can be obtained when 
the simulations involve approximately 10 grid nodes per minimal 
S-wavelength at the frequency /max of the source wavelet, even in 
the case of complex topographies with variable curvatures. 

4.5 Test 4: long-term stability 

The fourth test focuses on long-term stability jStacevl 1 19941 : 
iHestholni lioo i ). For this purpose, we consider a circular elastic 
domain with a radius of 150 m, surrounded by vacuum. The source 
S is located inside the circle, at (320 m, 200 m). This configuration 
is obviously not realistic, but it enlights the influence of the bound- 
ary on the numerical solution after many reflections, and especially 
on the possible excitation of numerical spurious modes leading to 
long-term instability. The mesh size is /i = 1 m. Time integration 
is performed during 10® time steps, with /max = 160 Hz. 

Figure[TOl-(a) shows a snapshot of at the final instant: no in- 
stability is observed, and the antisymmetry of is satisfied. Once 
the source is extincted (t > 2 tc), the mechanical energy E is theo- 
retically maintained. It can be written in terms of v and cr 



A 



2^i(A + ^) 



(20) 



At each time step, the integral in l |20| l is estimated by a basic trape- 
zoidal rule at the grid nodes inside Q. Figure[T0l(b) shows the time 
history of this mechanical energy so-obtained. It slightly decreases, 
due to the numerical diffusion of the scheme, which confirms that 
the method is stable. 



5 CONCLUSION 

Here we have presented a method of incorporating free boundaries 
into time-domain single-grid finite-difference schemes for elastic 
wave simulations. This method is based on fictitious values of the 
solution in the vacuum, which are used by the numerical integra- 
tion scheme near boundaries. These high-order fictitious values ac- 
curately describe both the boundary conditions and the geometrical 
features of the boundaries. The method is robust, involving negli- 
gible extra computational costs. 

Unlike the vacuum method, the quality of the numerical solu- 
tion thus obtained is almost independent of the angle between the 
free boundaries and the Cartesian meshing. Since the free bound- 
aries do not introduce any additional artefacts, one can use the same 
discretization as in homogeneous media. Typically, when a fourth- 
order ADER scheme is used on a propagation distance of 50 min- 
imal wavelengths, 10 grid nodes per minimal S-wavelength yield 
to a very good level of accuracy. With 5 grid nodes per minimal 
S-wavelength, the solution is less accurate but still acceptable. 

For the sake of simplicity, we have dealt here with academic 
cases, considering two-dimensional geometries, constant physical 
parameters, and simple elastic media. Let us examine briefly the 
generalization of our approach to more realistic configurations: 

(i) Extending the method to 3-D topographies a priori does not 
require new tools. The main challenge will concern the computa- 
tional efficiency of parallelization. A key point is that the determi- 
nation of each fictitious value is local, using numerical values only 
at neighboring grid nodes. Particular care will however be required 
for fictitious values near frontiers between computational subdo- 
mains, in order to minimize the exchanges of data. 

(ii) Near free boundaries, the domains of propagation are usu- 
ally smoothly heterogeneous. To generalize our method to continu- 
ously variable media, the main novelty expected concerns the high- 
order boundary conditions detailed in section 12.21 With variable 
matrices A and B indeed and k > 2, the procedure ^ will in- 
volve the following quantities, to be estimated numerically: 

gk-l Qk-l 

r— A, — B, Q = 0,...,fc-1. 

(iii) Realistic modeling of wave propagation requires to incor- 
porate attenuation. The only rheological viscoelastic models able to 
approximate constant quali ty factor over a frequency range are the 
generalized Maxwell body (Emmeric h & Korn[|l984h and the gen- 
eral ized Zener body tCarcione.,2001i) . These two equivalent mod- 
els ( IMoczo & Kristekl 20051) yield to additional unknowns called 
memory variables. In the time domain, the whole set of unknowns 
satisfies a linear hyperbolic system with source term 



= A-^U + B-^U -SU, 
at ox oz 



(21) 



where S' is a definite positive matrix. Compared with the elastic 
case ([TJ examined in the present paper, the main difference ex- 
pected concerns the time differentiation of the boundary condition 
l|2j. Indeed, equation ^ has to be modified accordingly to J21t . 
Similar modifications are als o foreseen in th e case of poroelasticity 
in the low-frequency range toaietal.L[T99"5) . where the evolution 
equations can be put in the form ( 121b . 
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